#!/usr/bin/perl -w
use strict;

use FindBin;
use lib ("$FindBin::Bin/..", "/net/cpp-group/Leo/bin", "$FindBin::Bin");
use db_parameters;
use ensembl_parameters;
use chromosomes;
use min_max;
use die_error;
use DBI;
use align_gaps;
use Data::Dump qw(dump);

print STDERR <<"HEADLINE";
888888888888888888888888888888888888888888888888888888888888888888888

    table_exon_sequences
	

	Description:
		Obtain the exonic DNA associated with each transcript
	Populates DB tables:
		prot_exon_sequences
		exon_sequences

888888888888888888888888888888888888888888888888888888888888888888888
HEADLINE

use constant GENEID =>0;
use constant CHRMSM =>1;
use constant START_ =>2;
use constant END___ =>3;
use constant STRAND =>4;
use constant CENTRE =>5;
use constant ORDER_ =>5;


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#	Determining exon structures of coding sequence


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888


#_________________________________________________________________________________________

#	open_chrm_file

#		retrieve genomic sequence

#_________________________________________________________________________________________
sub open_chrm_file($$)
{
	my ($taxon_data, $chrm_name) = @_;
	#
	#	return null if file not found
	#
	unless (open(DNA, "$taxon_data->[SEQ_DIR]/dna/chromosome.$chrm_name.fa"))
	{	warn
		"Error:\n\tCould not open '$taxon_data->[SEQ_DIR]/dna/chromosome.$chrm_name.fa'\n";
		return undef;
	}
	my $ignore = <DNA>;
	return (*DNA, length($ignore));
}



#_________________________________________________________________________________________

#	retrieve_chrm_sequence

#		retrieve genomic sequence

#_________________________________________________________________________________________
sub retrieve_chrm_sequence($$\$)
{
	my ($taxon_data, $chrm_name, $result) = @_;
	print STDERR " / Genomic DNA\n";
	#
	#	return null if file not found
	#
	unless (open FA, "$taxon_data->[SEQ_DIR]/dna/chromosome.$chrm_name.fa")
	{	warn
		"Error:\n\tCould not open '$taxon_data->[SEQ_DIR]/dna/chromosome.$chrm_name.fa'\n";
		return undef;
	}
	my $ignore = <FA>;
	my $seq = '';
	local $/ = undef;
	$$result = <FA>;
}


#_________________________________________________________________________________________

#	retrieve_exon_coordinates

#_________________________________________________________________________________________
sub retrieve_exon_coordinates($$$)
{
	my ($dbh, $taxon, $chrm) = @_;
	return $dbh->selectall_arrayref(<<"PL/SQLCMD");
	SELECT
		prot_id,
		prot_exons.start,
		prot_exons.finish,
		strand,
        exon_id
	FROM
		taxon.prot_exons NATURAL JOIN
        ens_id JOIN 
		taxon.gene_loci USING (gene_id)
	WHERE
		chromosome = '$chrm' and
        species = '$taxon'
	ORDER BY 
        gene_position,
		prot_id, 
		prot_exons.start;
PL/SQLCMD
}	


#_________________________________________________________________________________________

#	save_prot_coding_structure

#		retrieve exon coding structures for all proteins in a species

#_________________________________________________________________________________________

sub save_exon_sequences($$$)
{
	my ($dbh, $taxon, $taxon_data) = @_;
	print STDERR "\t$taxon Exons\n";

	#
	#	iterate chromosome at a time
	#
	my @chromosomes = get_chromosomes_with_sequences($dbh, $taxon, "$taxon_data->[SEQ_DIR]/dna");

	#DEBUGG
	#@chromosomes = (4);
	
	print STDERR "\t\t", scalar @chromosomes, " chromosomes for $taxon\n";

	my $verbose = scalar @chromosomes < 50;
	for my $chrm(@chromosomes)
	{
		print STDERR ("\t\tChromosome $chrm: Reading ") if $verbose;
		#
		#	exon coordinates ordered by protein
		#
		my $exon_coordinates = retrieve_exon_coordinates($dbh, $taxon, $chrm);
		unless ($verbose)
		{
			next unless scalar @$exon_coordinates;
			print STDERR ("\t\tChromosome $chrm: Reading ");
		}

		#
		#	chromosomal DNA sequence
		#
		#my $chromosome_seq;
		#retrieve_chrm_sequence($taxon_data, $chrm, $chromosome_seq);
		#next unless $chromosome_seq;
		my ($chromosome_FH,  $file_offset) = open_chrm_file($taxon_data, $chrm);
		next unless $chromosome_FH;

		
		print STDERR "(", 
						(scalar @$exon_coordinates), " exons / ";
		# for last exon of last protein
		push(@$exon_coordinates, ["LASTLAST"]) if @$exon_coordinates;


		my $cnt_seq = 0;
		my $curr_prot_id = '';

		# exon sequence making up current protein
		my %exon_ids;
		my @this_prot_exons_seq;
		my %prot_exons_seq;

		
		$dbh->do("COPY taxon.exon_sequences".
					"(exon_id, dna_sequence) FROM STDIN");

		# for each exon
		for my $line(@$exon_coordinates)
		{
			my ($new_prot_id,  $start, $finish, $strand, $exon_id) = @$line;

			if ($curr_prot_id ne $new_prot_id || $new_prot_id eq 'LASTLAST')
			{
				#
				# Start of new protein. Analyse exon structure of old protein
				#
				if ($curr_prot_id)
				{
					$cnt_seq++;
					print STDERR "." if ($cnt_seq % 10000 == 0);
					$prot_exons_seq{$curr_prot_id} = join ("",  @this_prot_exons_seq);
					@this_prot_exons_seq = ();
				}
				last if ($new_prot_id eq 'LASTLAST');
				$curr_prot_id = $new_prot_id;
			}

			# get DNA sequence of exon
			#my $exon_seq = substr($chromosome_seq,  $start - 1,  $finish - $start + 1);
			my $exon_seq;
			seek($chromosome_FH, $start - 1 + $file_offset, 0);
			read($chromosome_FH, $exon_seq,  $finish - $start + 1);
			$exon_seq = uc $exon_seq;
			die "$exon_id $start - 1 + $file_offset: ".($start - 1 + $file_offset). " -> ". ($finish - $start + 1) . "\n"
				unless length($exon_seq);


			#
			#	3' goes the other way: reverse and complement
			#
			if (!$strand)
			{
				$exon_seq =~ tr/ACGTacgt/TGCAtgca/ unless ($strand);
				$exon_seq = reverse $exon_seq;
				unshift(@this_prot_exons_seq, $exon_seq);
			}
			#
			#	Save 5'
			#
			else
			{
				push(@this_prot_exons_seq, $exon_seq);
			}

			# save exon sequence unless already saved
			unless (exists $exon_ids{$exon_id})
			{
				$dbh->pg_putline("$exon_id\t$exon_seq\n");
			}
			# save exon sequence
			++$exon_ids{$exon_id};
		}
		$dbh->pg_putline("\\.\n");
		$dbh->pg_endcopy;

		print STDERR "$cnt_seq transcripts)\n";

		$dbh->do("COPY taxon.prot_exon_sequences".
					"(prot_id, exonic_dna_sequence) FROM STDIN");
		$dbh->pg_putline("$_\t$prot_exons_seq{$_}\n") for (keys %prot_exons_seq);
		$dbh->pg_putline("\\.\n");
		$dbh->pg_endcopy;


		close $chromosome_FH;
	}
}







#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

# Main logic

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888



my $dbh = connect_to_panda();
use get_psql_time;
my $start_time = get_psql_time($dbh);

#
#	Get DB parameters
#
my %taxon_data = get_taxon_data;




#
#	taxon.prot_coding_exons
#
$dbh->do("TRUNCATE taxon.prot_exon_sequences");
$dbh->do("TRUNCATE taxon.exon_sequences");
$dbh->do("VACUUM FULL taxon.prot_exon_sequences;");
$dbh->do("VACUUM FULL taxon.exon_sequences;");








for my $taxon(keys %taxon_data)
{
	my (%prot_id_to_coding_exon_loci);

#	next unless $taxon eq 'monodelphis';
#	print STDERR "****************WARNING ONLY RUNNING FOR MONODELPHIS\n";
#    #
#    # Exons involved in coding per acc code
#    #
#	$dbh->do("SET WORK_MEM=1000000;");
#	$dbh->do("DELETE FROM taxon.prot_exon_sequences USING ens_id ".
#			 "WHERE prot_exon_sequences.prot_id = ens_id.prot_id AND species = '$taxon'");
#	$dbh->do("DELETE FROM taxon.exon_sequences USING ens_id, taxon.prot_exons ".
#			 "WHERE prot_exons.prot_id =ens_id.prot_id AND ".
#			 "exon_sequences.exon_id = prot_exons.exon_id AND species = '$taxon'");
	save_exon_sequences(	$dbh, 
								$taxon, 
								$taxon_data{$taxon});
}

$dbh->do("ANALYSE taxon.prot_exon_sequences;");
$dbh->do("ANALYSE taxon.exon_sequences;");


print STDERR "\tCompleted\n";


log_pipeline_stage($dbh, 'NULL', 21, 'table exon_sequences', $start_time);
$dbh->disconnect();

#remove zero size (empty error/status) files
#(-z $_) && unlink $_ for (glob("$dir_pipeline_errors/*"));

